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Computer simulations of realistic ion channel structures have always been challenging and 
a subject of rigorous study. Simulations based on continuum electrostatics have proven to 
be computationally cheap and reasonably accurate in predicting a channel’s behavior. In 
this paper we discuss the use of a device simulator, SILVACO, to build a solid-state model for 
KcsA channel and study its steady-state response. SILVACO is a well-established program, 
typically used by electrical engineers to simulate the process flow and electrical character- 


Keywords: istics of solid-state devices. By employing this simulation program, we have presented an 


Ion channel alternative computing platform for performing ion channel simulations, besides the known 


KcsA methods of writing codes in programming languages. With the ease of varying the differ- 
Poisson—Nernst-Planck ent parameters in the channel’s vestibule and the ability of incorporating surface charges, 
SILVACO we have shown the wide-ranging possibilities of using a device simulator for ion channel 
simulations. Our simulated results closely agree with the experimental data, validating our 

model. 
© 2006 Elsevier Ireland Ltd. All rights reserved. 
1. Introduction residues on the wall, the surrounding ions, and the charges 


Biological ion channels are transmembrane proteins that reg- 
ulate the ion flux through a cell, and thereby play an inte- 
gral role in cellular signaling mechanisms [1]. Modeling of ion 
channels has a long history, but the lack of detailed structural 
knowledge hampered the progress in this field. With advance- 
ments in molecular biology and X-ray diffraction techniques, 
a lotis now known regarding the structural details of ion chan- 
nels [2]. 

Ions in solution execute random Brownian motion, and 
if an electric field is applied, they acquire drift velocity. It is 
the interaction between the ions and the electric field in the 
channel that decides the salient features of an ion channel. 
The effective electric field is a complicated interplay from var- 
ious sources including the membrane potential, the charge 
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induced in the membrane wall by these ions [3]. 

Among the varied approaches available for ion channel 
modeling, the most accurate representation of the system is 
provided by molecular dynamics (MD) simulations, where all 
the atoms in a system are treated explicitly [4,5]. However, 
MD simulations are presently limited to simulation times of 
not more than a few nanoseconds, due to the small timesteps 
(~femtoseconds) required to resolve individual ion trajecto- 
ries. Any reliable estimates of steady-state channel currents 
require simulation time durations of the order of milliseconds, 
thus preventing an experimental validation of MD simulations 
[6]. 

An alternative method is to follow the motion of individ- 
ual ions using Brownian dynamics simulations [7,8]. Here, 
it is assumed that the protein structure is fixed and water 
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molecules are treated as a continuum. Collisions between the 
ions and the surrounding water molecules are mimicked by 
random fluctuating forces plus an average frictional force. A 
major advantage of Brownian dynamic approach is the ability 
to allow direct simulation of ion-ion interaction and calculate 
an ion’s interaction energy with respect to the system during 
each time step [7]. But this has always been an involved, time- 
consuming process [9]. 

Using a mean field approximation, continuum electrostat- 
ics provides a rather simplistic alternative to ion channel mod- 
eling [10]. It involves solving the Poisson—Nernst-Planck (PNP) 
equations for the charge distribution inside the channel, and 
has proven to be very computationally cheap. Even though 
the PNP method uses a major approximation of reducing the 
ion-ion interactions to an interaction between an ion and a 
mean field, its results have shown to be in good agreement 
with the experimental data [11,12]. 

In this paper we present a novel method of simulating 
ion channels based on continuum electrostatics. We have 
attempted to represent the structure of a KcsA channel as a 
solid-state device, with its intrinsic material properties similar 
to that in the realistic case [13,14]. Simulations are performed 
on a solid-state device simulator, SILVACO International, pro- 
viding us with the steady-state behavior of the KcsA channel. 
The simulation results agree well with the available experi- 
mental data. 


2. Methods 
2.1. Poisson—Nernst-Planck (PNP) theory 


The PNP approach to ion permeation in membrane channels 
has been studied in numerous works [10,11,15]. The PNP equa- 
tions are similar to the drift—-diffusion equations that describe 
the carrier dynamics in semiconductors [16]. The flux] of each 
ion species is described by the Nernst-Planck (NP) equation, 
which combines the diffusion due to a concentration gradient 
with the potential gradient [16,17]: 


zen 


=-D(|Vvn — V. ul 
j=-D(vn+ Ftv) (1 
where D is the diffusion coefficient, z the valency of the ion, 
e the electronic charge, k the Boltzmann constant, T the tem- 
perature, n the concentration of the ion, and @ is the potential. 
Also, the Poisson’s equation is written as [16,17]: 


eoV[e(r) Vo(r)] = —) zen — pex (2) 


ions 


where e(r) is the dielectric constant at an axial distance, r the 
sum over allions gives the total mobile charge, ¢9 the dielectric 
constant of free space, and pex represents all the external fixed 
and induced charges. It is to be noted that, realistically, all the 
quantities in Eqs. (1) and (2) are varying along the ion chan- 
nel axis [12,13]. This makes it notoriously difficult to obtain 
any analytical solution of the PNP equations, besides in some 
very specific cases [10]. However, it is possible to solve Eqs. 
(1) and (2) iteratively and obtain self-consistent solutions for 
the potential, concentrations and fluxes of the ions along the 


channel axis [17]. Even though this paper deals with 1D PNP 
computations, it has been shown that the PNP equations can 
be extended to 3D with a more realistic channel representa- 
tion and even closer fits to the experimental data [12]. 


2.2. The SILVACO simulator 


SILVACO is a solid-state process and device simulator that 
has the general capabilities for numerical, physically based, 
2D and 3D simulations of semiconductor devices [18]. It has 
two components: ATHENA, which is a 2D process simulation 
framework with software tools for modeling semiconductor 
fabrication processes, and ATLAS, which predicts the electri- 
cal behavior of semiconductor structures and provides insight 
into the internal physical mechanisms associated with device 
operation. ATLAS produces outputs in three forms: the run 
time output as a guide to the progress of the simulations 
running, log files storing summaries of the electrical output 
information, and the solution files storing the data related 
to the values of the solution variables within the device. The 
numerous carrier transport models, the availability of differ- 
ent material substrates, and the ability to account for other 
physical effects (e.g. surface charges, lattice heating, genera- 
tion and recombination) have added greater flexibility to this 
device simulator. 

Semiconductor device operation is modeled in ATLAS by a 
set of one to six coupled, non-linear partial-differential equa- 
tions (PDE) [18]. ATLAS provides numerical solutions of these 
equations by calculating the values of the unknowns on a 
mesh of points within the device. The original, continuous 
model (similar to the PNP equations) is converted to a dis- 
crete, non-linear algebraic system by an internal discretization 
procedure. The sets of PDEs, the mesh and the discretiza- 
tion procedure determine the non-linear algebraic problem 
to be solved. This is solved using an iterative procedure that 
refines successive estimates of the solution. Different solution 
procedures (e.g. Gummel, Newton, or block iteration) exhibit 
different behavior with respect to convergence, accuracy, effi- 
ciency and robustness. 


2.3. KcsA ion channel simulations with SILVACO 


Based on its structure revealed by X-ray crystallography [2,13], 
we have built a KcsA channel model as shown in Fig. 1. 
The potassium channel is modeled here as a transmembrane 
lumen, surrounded by protein walls and with cylindrical reser- 
voirs of potassium and chloride ions at two ends [4-13]. The 
channel extends from 10 to 50 A, with a narrow selectivity fil- 
ter of radius 1.5A and length 12A and a wider segment of 
length 23 A. The selectivity filter extends into the extracellular 
space, while the wider segment tapers towards the intracellu- 
lar space. Each reservoir at either end has a radius of 40 A and 
a width of around 8 A. 

To replicate the shape of the KcsA channel, the structure is 
made of two materials: one to mimic the water continuum that 
is conducting, and the other to mimic the non-conducting pro- 
tein walls [15,19]. The dielectric constants ¢ of the conducting 
and non-conducting mediums are set as 80 and 2, respectively. 
The entire device is constructed from 55 rectangular regions, 
which is the maximum limit for ATLAS [18]. Each conduct- 
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Fig. 1 - The model of a KcsA ion channel built in SILVACO. 
The channel extends from 10 to 50 A, with a narrow 
selectivity filter of radius 1.5 A and length 12 A and a wider 
segment of length 23 A. The selectivity filter extends into 
the extracellular space, while the wider segment tapers 
towards the intracellular space. 


ing region has the capability of being defined with its unique 
material properties. To this end, the carrier doping and their 
diffusion coefficients are specified in each conducting region 
as shown in Fig. 2. This allows us to emulate the fact that there 
are considerably fewer carriers in the vestibule of the channel 
compared to in the reservoirs. Interestingly, the selectivity fil- 
ter hardly has four to five ions at one instant [13], which is 
represented as an undoped region in Fig. 2. After the defini- 
tion of the device structure, the mesh points are specified to 
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Fig. 2 - The doping profile of the various conducting 
regions in the KcsA model. The extracellular baths have a 
symmetric concentration of 100 mM equivalently, while the 
selectivity filter is left undoped. The profile to chosen to 
represent the realistic scenario in the vestibule of the 
channel [13]. 


denote the locations where solutions are to be determined. 
The mesh is defined by a series of horizontal and vertical lines 
and the spacing between them. With a maximum limit of only 
9600 mesh nodes in ATLAS [18], it is important to use the mesh 
lines wisely and in locations which are crucial for the device 
behavior. Two electrodes, the source and drain, each of length 
20A are placed at either ends of the reservoirs for applying a 
potential and measuring the electrodiffusion characteristics 
of the device. 

Despite the fact that ion channels can be modeled as 
nanoscale electronic devices as discussed above, there are sev- 
eral issues to be addressed in their electrical simulations. The 
carriers in electronic devices are electrons and holes, while 
those in ion channels are charged ions bound by hydration 
shells [10]. The masses and radii of the semiconductor carriers 
are significantly small compared to those of ions in channels. 
Furthermore, the diffusion coefficient of electrons is around 
500cm?s~!, while that of K* ions is around 5 x 107-7 cm? s7!. 
The energy band structure and crystalline structures of solid- 
state materials is well known, whereas that of liquids is less 
known. On the brighter side, the drift-diffusion PNP theory still 
holds for both solid-state devices and ion channels, providing 
a simplistic, macroscopic study of the flow of electrons and 
ions alike [15-17]. With this basis, we have incorporated many 
of the characteristics of ion flow into our channel model. Even 
though it is difficult take into account the large mass and radii 
of ions, we have adjusted the diffusion coefficients of elec- 
trons in each conducting region to describe the slower motion 
of ions in a channel’s vestibule [19]. Also, the net carrier doping 
and the degree of ionization is specified in a manner to pre- 
dict decreased concentrations in the channel’s vestibule and 
to obtain realistic current-voltage plots, as shown in Fig. 2. 
For this purpose, we have used the fact that the carrier doping 
density (per cm?) is given by n= 10° NayC, where Nay is the 
Avogadro’s number and C is the concentration of the solution 
(mol/l). As will be demonstrated in the next section, with these 
modifications and appropriate mesh definition, it is possible 
to obtain realistic simulations results. 


3. Results 


The solid-state KcsA channel model was simulated using 
the PNP theory under various symmetric bath concentrations 
and potentials. The discretized PNP equations were solved by 
ATLAS using Gummel iteration and with each simulation run- 
time of <30s. 

Figs. 3 and 4 show the potential contours along the chan- 
nel at a bath concentration of 100 mM and voltages of 100 and 
200 mV, respectively. In Fig. 3, for example, a voltage potential 
of 100 mV was applied to the drain electrode in the reservoir 
at the right side, while the source electrode at the left side 
reservoir was kept at ground. The potential drops gradually 
along the channel, with most of the significant voltage drop 
occurring within the 12A length of narrow selectivity filter. 
Figs. 5 and 6 plot the electric field contours along the channel 
axis at 100 and 200 mV, respectively, and a 100 mM bath con- 
centration. The electric field variations can be observed along 
the channel, with noticeable changes at either mouths of the 
narrow selectivity filter. The entrance and exit of the channel 
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Fig. 3 - The SILVACO contour profile of the potential 
variation along the channel axis with 100 mV applied 
across the channel and extracellular bath concentrations of 
100 mM. 
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Fig. 5 - The SILVACO contour profile of the electric field 
variation along the channel axis with 100 mV applied 
across the channel and extracellular bath concentrations of 
100 mM. 


are made rounded, as any sharp boundaries cause problems 
while solving the Poisson’s equation and lead to unwanted 
peaks. 

Figs. 7 and 8 plot of the electric fields and potentials, respec- 
tively, along the channel for varying drain voltages and a bath 
concentration of 100mM. As observed, the simulation pro- 
gram takes into account an inherent built-in potential of the 
structure, which is around ~0.575V in this case. The origin 
of this potential is similar to that of Nernst Potential in ion 
channels as the SILVACO program considers the same Boltz- 
mann energy distribution of carriers [10]. In this case, the bath 
concentration is 6 x 10!9?cm~? (or 100mM equivalently) and 
the inside of the channel has an intrinsic concentration of 
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Fig. 4- The SILVACO contour profile of the potential 
variation along the channel axis with 200 mV applied 
across the channel and extracellular bath concentrations of 
100 mM. 


101° cm~3. Plugging this in the formula for the built-in poten- 
tial V=(kT/q) In(Cbath/Cchannel), WE get a value of ~0.575V. In 
Figs. 9 and 10, we plot the electric fields and potentials, respec- 
tively, along the channel for a fixed drain voltage of 100 mV 
and varying symmetric bath concentrations. As expected, the 
built-in potential increases with increasing extracellular con- 
centrations. 

In an attempt to check the validity of our model, we com- 
pare the current-voltage (IV) characteristics from our simula- 
tions with those of experiments [20]. Fig. 11 shows the corre- 
sponding IV plots of a KcsA channel at two symmetric bath 
concentrations of 250 and 500mM. Using reverse engineer- 
ing, the material properties of the solid-state KcsA structure 
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Fig. 6 - The SILVACO contour profile of the electric field 
variation along the channel axis with 200 mV applied 
across the channel and extracellular bath concentrations of 
100 mM. 
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Fig. 7 — Plots of electric field variation along the channel 


axis at varying voltages applied across the channel and 
extracellular bath concentrations of 100 mM. 


have been adjusted so as to mimic the realistic ion-channel 
structure and obtain a good fit between the IV curves. As seen 
from Fig. 11, the simulated IV curves closely agree with the 
experimental IV curves until 70 mV. Above this voltage, there 
is a slight tendency towards saturation in the experimental 
curves, whereas the simulated IV curves keep rising almost 
linearly. Nevertheless, this degree of agreement between our 
simulated results and the experimental IV curves is much bet- 
ter than those reported in the previous literature [9]. 

Next, we consider the effect of surface charges on the 
electrical properties of KcsA channels. These surface charges 
are known to influence the gating, conductance, and toxin- 
binding effects of KcsA ion-channels, and their actual physical 
locations are presumed to be in the selectivity filter [21]. Being 
negatively charged, these surface charges aid the permeation 
of cations and impede the flow of anions through the KscA 
channel. Here, we explored the effects of such charges on ion 
permeation by placing three positive charges in the selectiv- 
ity filter at the interface of water and protein walls (at 40, 44, 
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Fig. 9 — Plots of electric field variation along the channel 
axis at varying extracellular bath concentrations and an 
applied voltage of 100 mV across the channel. In each case, 
both the baths are symmetric, having the same 
concentration. 


and 48A in Fig. 2). It is worth mentioning that the surface 
charges in our simulations are positive (each Qsc = 10'4 cm~*) 
since the majority carriers used in our model are electrons. 
The SILVACO program is versatile in allowing us to incorpo- 
rate any fixed charges, interface traps, or even bulk traps at the 
semiconductor-insulator interface and observing their effects 
on the conduction of a device [18]. Fig. 12 plots the poten- 
tial variation along the channel with a bath concentration 
of 100mM and drain voltages of OV, 150mV and 250mvV. In 
comparison to its steadily increasing nature in Fig. 7 (with no 
charges), the axial potential variation in Fig. 12 (with surface 
charges) has a plateau inside the selectivity filter, with three 
small humps resulting from the charges at those locations. 
The simulated IV curves (with and without charges) are plot- 
ted in Fig. 13 for a 250mM bath concentration. Besides the 
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Fig. 10 — Plots of potential variation along the channel axis 
at varying extracellular bath concentrations and an applied 
voltage of 100 mV across the channel. In each case, both the 
baths are symmetric, having the same concentration. 
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Fig. 11 - The current-voltage characteristics of the KcsA 
channel model at two different symmetric bath 
concentrations (250 and 500 mM). Comparison is made 
between the SILVACO simulated plots (—, 250 mM; ---, 
500 mM) and the experimental data (A, KcsA: 250mM; ©, 
KcsA: 500 mM). 


general increase in conductance, there is an added asymme- 
try in the linearity of the IV curves (in positive and negative 
voltage ranges) by surface charges. As seen in Fig. 13, the 
current component due to the majority carriers (electrons) is 
enhanced much more compared to that due to the minority 
carriers (holes) owing to the attractive and repulsive forces, 
respectively, exerted by the surface charges. 


4. Discussion 


We have demonstrated the possibility of representing a KcsA 
channel as a solid-state device, and employing a device simu- 
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Fig. 12 — Plots of potential variation along the channel axis 
at varying voltages applied across the channel, in the 
presence of surface charges, and extracellular bath 
concentrations of 100 mM. Three pairs of charges (each 
Qsc = 1014 cm~?) are located inside the selectivity filter at 
the interface between the water and protein walls (at 40, 
44, and 48 A in Fig. 2). 
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Fig. 13 - The simulated current-voltage characteristics of 
the KcsA channel model in the presence of surface charges 
and at two different symmetric bath concentrations (250 
and 500 mM). Besides the increase in conductance, there is 
an enhancement in the majority carrier current and 
decrease in the minority carrier current. The nature and 
amount of surface charges is same as that in Fig. 12. 


lator, SILVACO, to simulate its electrical characteristics. With 
a goal to mimic the realistic and computationally feasible 
scenario inside an ion channel, we created a model with 
well-adjusted material parameters and employed SILVACO to 
obtain self-consistent solutions of the axial potential and ion 
fluxes. The simulated IV curves closely match the available 
experimental results, supporting the validity of the model. 

Besides giving insight into the general ion permeation 
mechanisms in a channel’s vestibule, the simulation results 
do open exciting opportunities to extend this work. Addi- 
tional features can be incorporated in the model to bridge the 
gap between the nature of a semiconductor-insulator and a 
water-—protein wall interface. In solid-state devices, the inter- 
face between a semiconductor and an overlying insulator is 
sharp and abrupt, while the dielectric boundary at a channel’s 
water-protein wall interface is not so distinctly defined. This 
dielectric boundary and its curvature is important in decid- 
ing an ion’s potential energy, which varies as 1/ewater in the 
baths and as 1/epotein Inside the channel [22]. Even though 
the precise value of éprotein is not crucial, it would better our 
model to put an additional protein region (€protein = 10) at the 
dielectric boundary to smoothen the sharp water-protein wall 
interface. Also, it is suspected that the dielectric constant of 
water decreases inside the channel, and recently Monte Carlo 
simulations have shown that lowering its dielectric constant 
in a channel’s vestibule alters a channel’s selectivity for dif- 
ferent ions [5]. Studies can be performed to investigate the 
effects of changing ewater on the channel permeation and its 
selectivity for monovalent or divalent ions. 

Although our model is based on macroscopic drift- 
diffusion computations, they do describe the experimental IV 
data surprisingly well. It is, however, tempting to be able to 
perform simulations at molecular and atomistic levels (Monte 
Carlo and Brownian dynamics) with the computational ease 
of PNP programs [19]. In this context, a common platform can 
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be built which has the capability of running molecular-level 
simulations for a short timeframe, inferring the ion-transport 
parameters, and feeding them to the SILVACO program to cal- 
culate the ion fluxes within reasonable time. Such a hybrid 
model would blend the accurate predictions of molecular-level 
simulations with the computationally cheap macroscopic flux 
calculations. 
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